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THE STATISTICAL ANATOMY OF OCEAN WAVE SPECTRA 


by 


Leon E. Borgman 


I. INTRODUCTION 


In the future, more and more ocean engineering design considerations 
will involve the wave energy spectrum. Typically, various secondary 
calculations will be made from the spectral and cross-spectral estimates 
at and between various space locations. The statistical reliability of 
the values obtained from the secondary calculations depend critically on 
the inherent statistical variability of the spectral estimates. 


The estimation of directional wave spectra is obtained by just such 
secondary calculations on the auto- and cross-spectral densities or 
corresponding finite Fourier transform coefficients for various wave 
properties measured at one or more space locations. The reliability of 
the directional spectrum depends both on the method of computation and on 
the intrinsic statistical variability of the Fourier coefficients or spec- 
tral estimates. Such characterizing quantities as the main direction of 
wave travel for a given wave frequency or some measure of the arc of 
directions from which waves of a given frequency are coming each have their 
own confidence intervals which ultimately relate back, through the method 
of calculation, to the spectral and Fourier coefficient variability. 


Over the years various theoretical probability relations have been 
derived which apply to linear waves (Pierson, 1955; Goodman, 1957; 
Blackman and Tukey, 1958). However, engineers are usually concerned with 
wave heights large enough to make linear assumptions questionable. Waves 
in hurricanes and other severe storms are prime examples of this situation. 
Yet it 1S just in such situations where probability confidence statements 
for the wave spectra or for derived secondary quantities are needed. 


In the following report, the statistical variations in wave energy 
spectral estimates for hurricane waves are examined empirically for 12 
separate intervals of wave record measured during Hurricane Carla (Septem- 
ber 1961). The measurements were made on a Chevron 0il Company platform 
in South Timbalier Block 63, Gulf of Mexico, in a 100-foot water depth. 
Hurricane waves were chosen for the analysis because they would illustrate, 
in exaggerated form, the effects of departures from linearity on the 
Statistical variability in spectral estimates. 


Various aspects of the study were reported at the 13th International 
Conference on Coastal Engineering in Vancouver, B.C., 10 to 14 July 1972, 
for one 20-minute record (Borgman, 1972). This study gives the analysis 
for the whole storm and develops certain implications and consequences of 
the empirical results. 


II. THE SPECTRAL COMPUTATION 


Two basic methods have been used in the past for computing the wave 
spectrum. The earlier one was based on the covariance function which was 
then numerically Fourier transformed and smoothed. That is, if Nn 
(n = 0,1,2..., N-1) is the water level elevation above mean water level, 
then the covariance function is: 


N-k-1 


Neuter 
Gh SS ear Ss "A"n+k ° @) 


(The quantity, N, is 4,096 for the analysis of Hurricane Carla.) Usually 


ey would be negligible for k larger than some value, say k,. Thus, 
the N numbers would be adequately summarized in the k,+1 values of Cy, 


Ordinarily, km is selected to be around one-tenth of N for most 
analysis of this type. The spectral density would be obtained from the 


numerical transform of the Cy: 


a Sie R mn 
p(fy) = At [& > DO De oF cos a Cs cos rn | , (2) 
q=1 
where 
fT (3) 
for r=0,1,3,..., k, (Blackman and Tukey, 1958). The quantity, At, is the 


timelag between successive measurements of n,(At = 0.2 second for the 
Hurricane Carla data). 


The computation of equation (1) entails a loss of information (N values 


replaced by k, values). This causes p(f,) in equation (2) to be a smoothed 


version of the true spectral density and distorts or eliminates features 
of the spectrum. The method of computation prevents the user from seeing 
aspects of the spectrum which may be important. 


The second method for computing the spectral density, which has come 
into wide favor during the last few years, is based on the application of 
the fast Fourier transform computing algorithm to the water level eleva- 
tions (Borgman, 1973). Complex-valued Fourier coefficients, A,, are 
obtained (or ™m—OF 25a N=) bye: 


N-1 N-1 
277 ; .. 2T7mn : 
A, = At ye ky COS coo iAt » Mn) SOD sa We aD (Ct) 
n=0 n=0 


Thus, Um is the cosine transform while V, is the sine transform of the 
water level elevations. The spectral lines are then computed by the 
formula, 


P(fm) = (Um + Vin)/(NAt) , (5) 
where 
m 
fn = WE (6) 
for m=0,1,2,3,..., N-l. These N spectral lines can be computed with 


great rapidity on a digital computer with the fast Fourier transform (FFT) 
procedure (Cooley and Tukey, 1965; Robinson, 1967). 


The frequency 
fyy = 1/2At (7) 
is called the Nyquist frequency. A symmetry relation 
Pen) = DCENem) (8) 


holds because of intrinsic mathematical properties of equations (4), (5), 
and (6). Hence, it is only necessary to specify PGES) LOT SO <msN/- Zi Gane, 
for frequencies between zero and the Nyquist frequency). The spectral 
density defined in equation (5) is defined by analogy to the conventions 
used by Blackman and Tukey (1958) to be a two-sided spectral density in 
which only the right-hand side is reported. That is, the total variance 


Nevers 
a Capt. 


of the water level elevations would be equal to 2 

The second procedure, involving the fast Fourier transform, was 
selected for the spectral computations for Hurrican Carla because it 
incurred less loss of information than the covariance procedure. The FFT 
method permits the inspection of all 2,048 spectral lines for frequencies 
up to the Nyquist frequency before the lines are averaged to yield a 


smooth estimate, p(f), of the spectral density. The covariance method 
gives only the smoothed spectral density. 


However, conclusions concerning statistical confidence for spectral 
estimates based on data computed with the FFT method are valid also for 
spectral estimates based on the covariance method. The two procedures 
for computing the spectral density give essentially the same result. 


III. ESSENTIAL EQUIVALENCE OF THE FFT AND COVARIANCE METHODS 
The covariance procedure gives spectral estimates which approximate a 


smoothing of the population or "true'' spectral density. The effective 
width of the smoothing is approximately 1/(2k,At), where k,At is the 


maximum lag used in the covariance estimate (i.e., k, is the maximum value 


of k used in computing equation (1)) and At is the time increment between 
water level measurements (Blackman and Tukey, 1958). (See App. A for the 
definition of effective width.) If "hamming'' smoothing, 


pC) = WAC) = CAC) = CAC aaa) . (9) 


is then made, the effective smoothing width for the estimates, DACEDE 
is changed to 1/(k,At) (Blackman and Tukey, 1958). 


The FFT method also yields estimates which are smoothed versions of 
the true underlying population spectral density. This can be seen from 
the following derivation. 


The true or population spectral density is defined as the integral 
Fourier transform of the theoretical covariance function. That is, 


Ge) = JailinGe)) MCE sy (10) 
where E[*] denotes the expectation operator and 


oo) . 
-127£T 


p(f) =| C(t) e dt (11) 


(Blackman and Tukey, 1958). Since C(t) is symmetric about T = 0 for 
stationary stochastic processes, it follows that equation (11) reduces to: 


foe) 


p(f) - f C(t). Cos, (Zee) Ge. (12) 


Equations (1) and (2) are obvious analogues to equations (10) and (12). 
It is not immediately obvious that the same thing can be said for equa- 
tions (4) and (5), although it is true there also. To see this, let 


N-1 


a it 
Cy =. N Le Nn Nn+k ’ (13) 
n=0 


where N = 4,096 in the context of the Hurricane Carla records and the n 
for n>N and n<O required for the computation of equation (13) are 
defined by periodicity as: 


= Thay. (14) 


Equation (13) will give almost the same estimate of the covariance function 
as equation (1) for O<k<N/2. There is, of course, a little distortion 


related to the effects of equation (14). However, this will be over- 
whelmed by the averaging against the other "within sequence" lag products 
provided water level observations separated by more than k, time incre- 


ments are essentially independent of each other and k, is small. 


The periodicity introduced in equation (14) causes a corresponding 
periodicity in CL. 


CMC euCyay. for y.N/2<k<N (15) 


An approximation to equation (12) would be 
N/2 
Pea) tye NC Hicos) (27kkAe)y ar! (16) 


k=t +] 


With the introduction of equation (15) and the definition of f, (equa- 
tion 6), this becomes: 


N-1 N-1 

nw n~ m an 

p(fm) = At » Ck cos (2nsRe at) = At ye Ge cosm(2imk/N)) ven) l7) 
k=0 k=0 


/ 


The last equation is algebraically equivalent to equation (5) although 
it appears somewhat different. The equivalence can be seen after equation 


(13) is substituted for Cy, and the transform is shifted to exponential 


form. Thus, 


Z| 


N-1 N-1 

fs 1 a 

Ge) = We », ( Ny Mea exp (-i21mk/N) 
k=0 n=0 


N-1 N-1 
= ot yy In, exp (+4 2nmn/) [oe exp (-i2mm(ntk)/N] 
k=0 n=0 
N-1 Neal 
= = In, exp c2nan/X) i n exp (-idnmk/n| (8) 
n=0 k=0 


after the introduction of the periodicity assumed in equation (14). It 
follows that 


Gen) = (ny Aa) /NiNe S (UE WE)/NieY (19) 


where Am denotes the complex conjugate of Am. This is the desired con- 
clusion and equations (5) and (17) have been shown to be equivalent. 


Returning to the question of the amount of spectral smoothing involved 


in the estimate, DiGawe define the finite Dirac comb (Blackman and Tukey, 
1958) as: 


c-1 

Se At ‘ 

Vo (t3at) = 5 6 (t + ¢ At) + at De S(t - kat) + 5 6(t - ¢ At) , 
k=-c+l1 (20) 

where 6(x) is the Dirac function which has the properties: 
if SCgeb<e = i (21) 

SG) sO | if x#0 (22) 
and for any bounded function g(x): 

f ste - etdax = afc) - (23) 


Now suppose the sampling questions involved C, are ignored and C(kAt) 
1s inserted in) placenot ) Cy. an) equation) i@o)e- Dheny, 


foe} 


p(Eq) =i Ga) Yayo Cents) Ge Nake (24) 


The right-hand side of this equation is the Fourier transform of 
iG) Vn/2 (CegAe)l Thus, by Fourier transform theory p(f,) is the convo- 


lution of the transforms of the separate functions or 


eS le (25) 


=-0O 


P(fm) ~ p(£m) * 
k 


In the above p(f,) is the transform of C(T) by equation (11) and 
© k J ! 
He Qo (F - ie) is the transform of Vn/2 (t;At) (Blackman and Tukey, 


1958). The function Qo(f) is (Blackman and Tukey, with Ty = NAt/2): 


_ sin(Tf£N At) 
Zo = ——s (26) 


Hence, 


2 


B(f_) ; 


2, i Qh Eni x - y)p(y)dy . (27) 


This shows that DiC) to the extent that Cy behaves like C(kAt), is 
the result of smoothing with the Qjo(f) function and then summing or 
aliasing as indicated by the summation. 


The function Qo(f) has an effective width of (1/NAt). Hence, 9 (Ge) 


represents a smoothing approximately over a frequency interval of (1/NAt). 
This is the spacing between the f, in equation (6). It follows that 
2ach FFT spectral estimate represents approximately a smoothing of p(f) 


WereheIncernval t= A @y/2NAE): 


It should be noted that problems related to the side lobes of Qo(f), 
the aliasing, and the sampling variability have been ignored in the above 
jiscussion. These effects will now be discussed briefly. The distortion 
due to aliasing can be mitigated by choosing the Nyquist frequency suf- 
ficiently large. The sampling variability shows up in the statistical 
fluctuations of the estimates and can be examined from that perspective. 
If the true spectra are relatively smooth and linearly changing, the side 
lobes of Qg(f) will compensate for each other somewhat. However, if the 
true spectra are not smooth, there will be unavoidable leakage of large 
spectral departures or spikes into neighboring spectral estimates. 


IV. FFT AVERAGING NEEDED TO GIVE SPECTRA EQUIVALENT TO THOSE 
FROM THE COVARIANCE METHOD 


From the previous section, the covariance spectral estimates (after 
hamming) represent a smoothing over the frequency interval of width, 
(1/ky At), while the FFT spectral estimates involve a smoothing over a 
frequency interval of width, (1/N At). Hence, the number of FFT spectral 
lines that should be averaged together to yield an estimate with approxi- 
mately the same smoothing as the covariance spectral estimates is: 


: (1/km At) 


number = 
(1/N At) 


Psi (28) 


V. SUMMARY OF COMPUTATIONAL FORMULAS USED FOR THE HURRICANE CARLA 
WAVE SPECTRA 


Twelve pieces of record at various times in the storm (Table 1) were 
chosen for analysis. Each piece consisted of N = 4,096 values of water 
level elevations taken at a time increment of At = 0.2. Thus, each of 
these records was NAt = 13.65 minutes long. The frequency scale incre- 
ment would be Af = 1/NAt = 0.00122 sec.~!. The Nyquist frequency for 
this choice is fy = 1/2At = 1/0.4 = 2.5 sec.~} which corresponds to a 


period of 0.4 second. Energy in waves with periods smaller than 0.4 second 
would thus be aliased into lower frequencies. An examination of the final 
spectra produced little evidence of much energy past this selected Nyquist 
frequency. 


Table 1. Hurricane Carla data analyzed. 


Date Data code 


8 Sept. 1961 


9 Sept. 1961 


10 Sept. 1961 


The 4,096 water level elevations were transformed by the fast Fourier 
transform algorithm to yield Fourier coefficients, for O<m<N: 


4,095 
A, = At So. Tn eo 12mmn/N (29) 
n=0 
The spectral lines were then computed from 


DCE) S |Anle/ CG1uOs2). 5 (30) 


when | Am | denotes the complex modulus of A, . The spectral density was 
estimated by a moving average of the spectral lines: 


B(fm) = Ei) B(f-5)/ ae (31) 


where 


W5 = exp (Gas /ats) ae (32) 
These weights are, thus, Gaussian smoothers with a standard deviation of 

o = 3Af = 0.00366 sec.~!. The numerical values of the Wj are given in 
Table 2. In the vicinity of zero frequency, the subscript j in equation 
(31) was summed over the possible values and the divisor normed the weights 
appropriately. At m= 6, for example, j was summed over -6<j<13 (the 
left tail of the moving average was truncated off). we 


Table 2. Weights used in the moving average 
of the spectral lines to produce 
the estimates of the spectral density. 


au 
= 
Cat, 


0) 
1 
2 
3 
4 
5 
6 
ii 
8 


SOLO OOO Oo Oroiore or 
(>) 
lon 
ea) 
—N 


Note: w_; = w:. 


The effective width of the Gaussian smoother in equation (32) is V2T oO 
where oO = 3Af, if lags are measured on the frequency scale and o =3, 
if lags are measured relative to the number of spectral lines. The effec- 
tive width in terms of number of spectral lines encompassed is thus 


3 V2n = 7.52, or rounding to the nearest integer, 8 spectral lines. By 
equation (28), the FFT spectral density estimates so produced would be 
comparable to covariance spectral density estimates with maximum covariance 
lag derived from: 


N 
se 8} (33) 
kn 
or 
4,096 
ky = Se = S12 (34) 


VI. DISCUSSION OF HURRICANE CARLA WAVE SPECTRA 


The spectral lines (shown as dots) and the spectral density estimates 
averaged from the lines (shown as a line or a string of pluses) are given 
in Figures 1 through 12. 


Close examination yields two very interesting facts. First, the peak 
values of the spectral densities are related closely to one or two spec- 
tral lines. Only in record 6883 is the peak related to four exceptionally 
large lines. In the other cases it is always one or two. Second, these 
exceptionally large spectral lines do not seem to persist. The two 
members of record pairs (6881-1 and 6881-2; 6886-1 and 6886-2) are sepa- 
rated from each other by only 20 minutes. Yet in both cases, exception- 
ally large spectral lines are present in one member of the pair but not 
in) the others. 


A general examination of the spectral lines versus the spectral 
density estimates cannot help but develop a sense of healthy skepticism 
concerning the general reality of the fine structure in the spectral 
density. JAlson ite ais tellit thats the texceptionaliliy larce spectralaglamiess 
often twice as large as the nearest other line value, must belong to 
another population from the rest of the lines. Perhaps some sort of 
resonant phenomenon is creating a main wave train with the rest of the 
lines functioning as superimposed noise. 


VII. STATISTICAL VARIABILITY OF THE SPECTRAL LINES 


The spectral density was subtracted from each spectral line for 
6<m<305 to provide 300 residual values, R, , 


Rn = P(t) - B (£m) : (35) 


A positive and negative standard deviation for the residuals were computed 
with the weights introduced in Table 2. Let J, be the values of the 
index m for which the residuals are positive while J_ are the values 
of m for which the residuals are negative. The positive and negative 


variance of the residuals are defined as: 


2) ye -R2 . : 
OF eet 2 wees) / » wi (36) 
JES 5, jeJ, 
625 a w5Ra-3 / ee (37) 
Ain aM J m-) 4 J 
jeT_ jeJ_ 
where Mev emeansm bDelongs Comune Seta Oteuus wisi O) will be a moving 
+,m 


average estimate of the root-mean-square (rms) of the positive residuals 
in the vicinity of the frequency f,. A similar statement relative to 


O_ ,, and the negative residuals will also hold. 
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Figure 1. Hurricane Carla spectrum record number 6877, 
0600 hours, 8 September 1961. 
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Figure 2. Hurricane Carla data number 6878, 1200 hours, 
8 September 1961. Various smoothings of spectral data. 
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Figure 4. Hurricane Carla spectrum record number 6880, 
0000 hours, 9 September 1961. 
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Figure 6. Hurricane Carla spectrum record number 6881-2, 
0620 hours, 9 September 1961. 
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Figure 7. Hurricane Carla spectrum record number 6882, 
1200 hours, 9 September 1961. 
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Figure 8. Hurricane Carla spectrum record number 6883, 
1500 hours, 9 September 1961. 
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Figure 9. Hurricane Carla spectrum record number 6884, 
1800 hours, 9 September 1961. 
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Figure 10. Hurricane Carla spectrum record number 6885, 
2100 hours, 9 September 1961. 
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Figure 11. Hurricane Carla spectrum record number 6886-1, 
0000 hours, 10 September 1961. 
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Figure 12. Hurricane Carla spectrum record number 6886-2, 
0020 hours, 10 September 1961. 
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Symmetrically normed residuals (SNR) are then defined as: 


a SO, m > ae I > 


SNR = (38) 
Lain rer aii iea (allen 


Thus, the symmetrically normed residuals are ratios of the residuals to 
the rms positive or negative residuals, as the case may be. 


Originally, the analysis was made with rms residual, disregarding 
whether the residuals were positive or negative. This, however, led to 
ridiculous results particularly in the subsequent use in simulation. The 
use of different norming divisors for negative and positive residuals 
avoided these peculiarities completely. 


The residuals and the corresponding SNR's are plotted in Figures 13 
through 24. The symmetrical normalization seems to very adequately pro- 
duce a uniform cloud of points. There does not seem to be any tendency 
for the SNR's to be systematically large or small at any particular 
frequency. Generally about 60 percent of the SNR's fall below zero. 


As a check against sequential dependence among the SNR's, neighboring 
pairs of SNR's were plotted on a scatter diagram in two-dimensional space. 
The first member of the pair of SNR's was the x-coordinate while the 
second member was the y-coordinate for the plotted point. Any tendency fo 
big values to follow big values (or the reverse) would show up as a clus- 
tering tendency on such a plot. Complete independence, on the other hand, 
would show up as a uniform cloud of points. 


The plots of this type for the 12 Hurricane Carla data sets are given 
in Figures 25 through 36. The point scatter is really quite uniform for 
all of the records with no obvious dependencies showing up. However, it 
should be emphasized that this type of examination only reveals overall 
average dependencies. There may be dependencies between neighboring SNR's 
at certain frequencies which are counteracted by opposing dependencies 
at other frequencies. However, the earlier graphs (Figs. 13 through 24) 
would show any strong dependencies tied to frequencies if they were 
present. 


Everything considered, the analyses strongly support the conclusion 
that the spectral fluctuations have been successfully decomposed into a 
smoothed spectrum plus a constant (either O, or O_) times independent 
random noise. The smoothed spectrum and the constants are frequency 
dependent. The noise apparently does not depend on frequency. In symbols, 
this decomposition can be written: 


B(fm) = B(Em) + ¢ ° (SNR) , (39) 
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Figure 14. Hurricane Carla data number 6878, 1200 hours, 
8 September 1961. Residual analysis, Gaussian 
smoothing on spectral lines. 
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Residual analysis, Hurricane Carla data number 6879. 
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analysis, Hurricane Carla data number 6880. 
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Hurricane Carla data number 6881-1. 
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Figure 18. Residual analysis, Hurricane Carla data number 6881-2. 
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Residual analysis, Hurricane Carla data number 6883. 
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Hurricane Carla data number 6884. 
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Residual analysis, Hurricane Carla data number 6885. 
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Residual analysis, Hurricane Carla data number 6886-1. 
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Residual analysis, Hurricane Carla data number 6886-2. 
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Scatter diagram: Sequential pairs of symmetrically 
normed residuals, data number 6877. 
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normed residuals, data number 6879. 
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Scatter diagram: Sequential pairs of symmetrically 
normed residuals, data number 6880. 
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Scatter diagram: Sequential pairs of symmetrically 
normed residuals, data number 6881-1. 
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Scatter diagram: Sequential pairs of symmetrically 
normed residuals, data number 6884. 
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Scatter diagram: Sequential pairs of symmetrically 
normed residuals, data number 6885. 
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Scatter diagram: Sequential pairs of symmetrically 
normed residuals, data number 6886-1. 
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where 
oO, 5 abe GR S @ 


Cir { , (40) 
o if SNR < 0 


and SNR denotes the SNR regarded as a random variable. 
VIII. THE EMPIRICAL PROBABILITY LAW FOR THE SYMMETRICALLY NORMED RESIDUAL 
The cumulative distribution function for the SNR is defined as: 


FoNR (w) = P [SNR < w] , (41) 


where P[°] is the probability of the event specified within the brackets. 
This distribution function can be estimated from the 300 values of the 
SNR's for each record of Hurricane Carla. Let (SNR), be SNR's ranked 


in order of increasing size: 
(SNR) 1 < (SNR), < (SNR), Oe < (SNR) 369 . (42) 


A statistically reasonable estimate of Fonr (Wy) for 


We = (SNR), (43) 
is 
E (wap piso (44) 
SNR* k 301 


(Gumbel, 1954). Thus, a graph of Pon (wy) versus w, for k = 1,2,3,°**,300 


gives the distribution function estimate. The graphs are shown in Figures 
37 through 48. 


The corresponding probability densities may be obtained by differen- 
tiating the distribution function numerically. For the present study, 
this was done by selecting a band on the SNR axis which is 0.5 unit wide 
and fitting a least square line to all the ranked points lying within the 
band. The slope of the line is the probability density estimate assigned 
to the midpoint SNR value for the band. This was repeated for all 300 
possible midpoints on the SNR axis. The resulting probability densities 
are given in Figures 49 through 60. 


IX. EMPIRICAL PROBABILITY INTERVALS FOR p (fm) BY SIMULATION 


Suppose a new spectral density estimate, io 2.) is developed by 


simulation from equation (39) by the following procedure. A random num- 
ber, uniformly distributed on the interval, (0,1), is generated in the 
digital computer. One of the pieces of record is selected for study and 
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Figure 37. Cumulative distribution function, Hurricane Carla data 
number 6877, 0600 hours, 8 September 1961. 
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Figure 38. Cumulative distribution function, Hurricane Carla data 
number 6878. 
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Figure 39. Cumulative distribution function, Hurricane Carla data 
number 6879, 1800 hours, 8 September 1961. 
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Figure 40. Cumulative distribution function, Hurricane Carla data 
number 6880, 0000 hours, 9 September 1961. 
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Figure 41. Cumulative distribution function, Hurricane Carla data 


number 6881-1, 0600 hours, 9 September 1961. 
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Figure 42. Cumulative distribution function, Hurricane Carla data 
number 6881-2, 0620 hours, 9 September 1961. 
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Figure 43. Cumulative distribution function, Hurricane Carla data 
number 6882, 1200 hours, 9 September 1961. 
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Figure 44. Cumulative distribution function, Hurricane Carla data 
number 6883, 1500 hours, 9 September 1961. 
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Figure 45. Cumulative distribution function, Hurricane Carla data 
number 6884, 1800 hours, 9 September 1961. 
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Figure 46. Cumulative distribution function, Hurricane Carla data 
number 6885, 2100 hours, 9 September 1961. 
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Figure 47. Cumulative distribution function, Hurricane Carla data 
number 6886-1, 0000 hours, 10 September 1961. 
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Figure 48. Cumulative distribution function, Hurricane Carla data 
number 6886-2, 0020 hours, 10 September 1961. 
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Figure 49. Probability density, Hurricane Carla data number 6877, 
0600 hours, 8 September 1961. 
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Figure 50. Density function, Hurricane Carla data number 6878. 
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Figure 51. Probability density, Hurricane Carla data number 6879, 
1800 hours, 8 September 1961. 
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Figure 52. Probability density, Hurricane Carla data number 6880, 
0000 hours, 9 September 1961. 
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Figure 53. Probability density, Hurricane Carla data number 6881-1, 
0600 hours, 9 September 1961. 
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Figure 54. Probability density, Hurricane Carla data number 6881-2, 
0620 hours, 9 September 1961. 
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Figure 55. Probability density, Hurricane Carla data number 6882, 
1200 hours, 9 September 1961. 
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Figure 56. Probability density, Hurricane Carla data number 6883, 
1500 hours, 9 September 1961. 
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Figure 57. Probability density, Hurricane Carla data number 6884, 
1800 hours, 9 September 1961. 
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Figure 58. Probability density, Hurricane Carla data number 6885, 
2100 hours, 9 September 1961. 
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Figure 59. Probability density, Hurricane Carla data number 6886-1, 
0000 hours, 10 September 1961. 
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Figure 60. Probability density, Hurricane Carla data number 6886-2, 
0020 hours, 10 September 1961. 
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the corresponding cumulative distribution function figure is picked from 
Figures 37 through 48. The coordinate value on the vertical axis of the 
distribution function which equals the uniform random number is located. 
Reading horizontally from this coordinate value to the empirical distrib- 
ution curve and then down vertically to the SNR axis yields a random SNR 
value. This procedure is illustrated in Figure 37 by the dotted line. 
The uniform random number is 0.532. The corresponding random SNR value 
is -0.16. 


The distribution function for the SNR values obtained by this proce- 
dure will be identical to the graphed empirical distribution function 
Fonr(w). This follows from the following argument. The SNR, so developed 


are less than or equal to w if, and only if, the uniform random number 
is less than or equal to Foyp(w). This is true because the two numbers 


are tied together via the graphed curve. Hence, 


(45) 


> 


P [random SNR < w] = P [U < Fonp(w)] 


where U denotes the uniform random number. But by definition, the dis- 
tribution function for a uniform random number is: 


Ey Cy eek LU Se Se (46) 


Hence, returning to equation (45), 


P [random SNR < w] -Ip (U 


[A 


The above procedure is repeated for 300 independent uniform random 
numbers to obtain 300 random SNR values. These 300 SNR values are just 
as likely to have happened as the originally occurring values, provided 
the decomposition in equation (39) is accepted as valid and provided the 
independence assumption truly holds. 


e Hence, equation (39) can be used with the 300 SNR values and the 
DCE Ne O+,m and O_ im frequency functions to create a new set of 
spectral lines, pace) These spectral lines might just as well occurred 


as the original set if the random spectral fluctuations had accidentally 
gone that way. 


Finally, the 300 simulated spectral lines, Da Geen) are smoothed 
according to equation (31) to produce new simulated spectral densities, 


B* (fm) - 


The above procedure in its entirety was repeated 900 times for each 
of the 12 pieces of Hurricane Carla data. Thus, 900 statistically equiva- 
lent spectral densities were generated by simulation for each hurricane 
data record. 
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n 
a 


How much do the p*(f,) as a group differ from the spectral density, 
D(E_) > used in the simulation? The answer to this question was developed 
for the 10 frequencies, 0.070, 0.072, 0.074, 0.077, 0.083, 0.088, 0.101, 
0.132, 0.168, and 0.243 sec.~!,. The 900 values of p* (fm) were ranked for 


each frequency and the two values with ranks 45 and 855 were selected as 
estimates of the 5th and 95th percentiles for p*(f,). These percentile 


estimates are plotted versus frequency as dots in Figures 61 through 72. 


X. COMPARISON WITH CHI-SQUARED PROBABILITY INTERVALS FOR A (£m) 


If the sea surface is Gaussian, the spectral density will follow a 
probability law closely related to a chi-squared random variable with 16 
degrees of freedom (for the Hurricane Carla estimates) (Borgman, 1972). 
Symbolically, i 

nS a) |e 
mew Tie ee 


where p(f,,) denotes the true or population spectral density. Thus, if 
denote the 5th and 95th percentiles for a chi- 


2 2 
SG 0.05 8) 2 Sa20.05 
squared random variable with 16 degrees of freedom, then: 


2 2 
x Ge.) x PGE) 
; | 16,0.05 ae ait) < 16,0.95 m 
16 16 


= 0.90 5 G49) 


. D 2 . 
The interval levanes DG) ye 16 , X16 0.95 p(£) i] 16 ) thus provides 


a 90 percent probability interval for CGE DE However, p(fm) is not 
known. As an approximation, f(f,,) may be substituted for p(f,). An 
analogous approximation was made in the simulations. The resulting upper 
and lower limits are plotted as asterisks in Figures 61 through 72 versus 
each of the selected frequencies. The spectral density values p(fm) are 


shown in the figures as pluses. 


As was noted in an earlier paper based on an analysis of part of the 
data (Borgman, 1972), the chi-squared probability intervals do not differ 
excessively from the simulated probability intervals. They are both about 
the same, although there are substantial variations from record to record. 
The upper bounds show appreciably more scatter than do the lower bounds. 
Another comparison of the two kinds of probability intervals is given in 
Figure 73. The two ratios, 


simulation upper bound / (£m) 


and, 


simulation lower bound / pice) , 
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Figure 64. Probability intervals for p, Hurricane Carla data 
number 6880. 
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Figure 66. Probability intervals for A, Hurricane Carla data 
number 6881-2. 
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Figure 67. Probability intervals for p, Hurricane Carla data 
number 6882. 
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Figure 68. Probability intervals for 6, Hurricane Carla data 
number 6883. 
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number 6886-2. 
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are plotted versus frequency for all 12 data records. The average over 
the 10 frequencies for each record is plotted on the right of the graph 
above the frequency value, 4.0 sec.~!. For the chi-squared distribution 
the corresponding values are theoretically, 


95th percentile 6 (£m) / p(fm) = ‘Gis 0 95 / 16 = 1.64 (50) 
Sth percentile p(f,) / P(f_) = X26 5 af) i648 0.80 2 Giy 


These values are shown as dashlines in the figure. 


On the average the two types of probability intervals agree fairly 
well. 


XI. COMPARISON OF THE EMPIRICAL DISTRIBUTION FUNCTION OF THE SYMMETRI- 
CALLY NORMED RESIDUALS WITH THE CHI-SQUARED VERSION 


If the chi-squared distribution is reasonably valid for (fm) » it 


may also hold for the spectral lines, at least approximately. The theore- 
tical relation to be checked for validity is: 


SC / DD 296 0 2 (52) 
The chi-squared analogy to o? in equation (38) is: 


E | | Stem) = p (fn) } ; B(fm) > p (én) | 


theoretical of 


P(f,) 2 pig) 
= pty) & [{ sey - 1} Tt 
eh ie hae 
By exact analogy: 
theoretical o2 = ptm) | {x3 i at | x3 < 2 | ; (54) 


Since the probability density for Ne random variable is: 


ea W/2 
5 A acon ny 2\\(0) 
2 (w) = { i (55) 
2 0 , for w <0 


Us) 


and, 


vu 
ra | 

>< 

Nod 

Vv 

i) 
a 

i} 


i -w/2 
f —=— dw = Sa (56) 
2 


iT] 
a 

i} 
i) 


p | <2 | ane (57) 


it follows that the conditional densities needed to evaluate equations 
(53) and (54) are: 


eo W/2 
5. one iy & 
2 oom a 
f 2| 255 (Ww) = (58) 
X51 X57 
D2 
0 » otherwise 
oe W/2 
Si » for O<w<2 
P Zier =) 
Nivea, (W)) = (59) 
x2] x32 
0 , otherwise 


Thus, the quantities in equations (53) and (54) are given by: 


co 


2 f -w/2 
theoretical On ae (fm) if (w - D2 Ses dw 


Age 2 ze 
(60) 
2 
ef 279-62.) = {1.414 p(n) } 
2 
2 -w/2 
a 
theoretical o% = a / (we=2)) 2 e dw 
mt adices .) 2 
0 
(61) 


oil 2 
Y (+s) P?(fq) = { 0.6465 p(n} 


l-e 
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Finally, the chi-squared analogy to equation (38) will be: 


D(fm) - fy 
theoretical SNR = P(fm) = p(fm) . (62) 


cp (fy) 


where 


We ANA” esse GSGee)) eS tales) 


QO 
tl 


(63) 
0.6465 , bts pice) <p) 
In terms of the 1G random variable, equation (62) becomes: 
a B (fm) 
theoretical SNR = c i | Bem - i 
p (fm) 
2 
= -1 Let 
ee) Wye 2 (64) 
Ae 2 


Hence, the theoretical distribution function for the SNR variables 
(consistent with the chi-squared assumption) is: 


: i ‘dan fan 
theoretical FoNR (w) = P We {x3 - 2} <w 


Hl 
ae) 
(SS) 
>< 
NO NO 
[A 
NO 
fe) 
= 
+ 
No 
eed 


= F 9 (2cw + 2) = 1 -exp(-cw - 1) 
X9 


l= exp(-1.414 w- 1) ; df w > 0 
= . (65) 
1 - exp(-0.6465 w- 1), if w< 0 


The theoretical F our) is graphed in Figures 37 through 48 as a solid 
line. 
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By inspection, it can be seen that distribution function derived from 
the chi-squared probability law is reasonably close to the empirical 
distribution functions although there are systematic differences. In 
general, the empirical curves tend to lie below the theoretical curve for 
most argument values. The Kolmogorov confidence interval for the dif- 
ference between the true distribution function and the empirical one is 
drawn in on record 6878. The chi-squared curve exceeds the upper boundary 
in the midranges but the two distribution curves are in fair agreement in 
the vicinity of the 5th and 95th percentiles. This is probably why the 
probability intervals agree fairly well. It should be noted that record 
6878 is one of the more extreme cases and most of the other records 
attain closer agreement between the two curves. 


A numerical error was made in the earlier report (Borgman, 1972, Figs. 
7 and 8) relative to the chi-squared related distribution function and 
probability densities. The theoretical curves in those figures should be 
ignored. 


XII. THE OUTLIER SPECTRAL LINES 


In 7 of the 12 records, one or more spectral lines loom high above 
the others. These are ordinarily associated with the largest spectral 
density values and determine where the peak of the spectral density will 
occur in most cases. The spectral density decreases appreciably if these 
extreme or outlier spectral lines are deleted from the averaging process 
in the density determination. 


A list of all the outliers is given in Table 3. Spectral lines which 
exceeded 500 are shown and one value which went to 477 is included. The 
spectral density value at that same frequency is tabled as well as the 
spectral density which results if the outlier spectral lines for that 
record were all deleted from the. averaging. The spectral density without 
the lines is usually about 60 percent of what the density is with the lines 
included in the averaging. 


Two other measures of the "extremeness" of the lines are listed in 
Table 3. The first is the ratio of the spectral line value to the density 
computed with the outliers deleted. If chi-squared theory holds, this 


ratio should behave like a y2 / 2 random variable. The 99.5 percentile 
X9 


for a x5 / 2 random variable is 5.3. Six of the 14 outliers listed 


exceed 5.5 in) values “However, itis ditticullt to antexrprets chisem Ones 
examining the larger members of 300 lines. Hence, extremal statistical 
theory needs to be introduced. Straightforward application of the theory 
of extremes would say that the probability that the largest such ratio 

in a record would be less than 5.3 is: 


O98)! = 0.22, (66) 


providing the chi-squared interrelation and the independence assumptions 
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Table 3. Characteristics of outlier spectral lines. 
Outlier 
Record No. Freq. line 
(a) 
6878 0.0732 477 148 96 
0.0842 525 WILY 54 
6879 0.0757 620 141 67 
6881-1 0.1025 652 231 166 
6882 0.0720 557 299 203 
0.0745 660 338 251 
6883 0.0696 586 i2 102 
0.0781 741 248 132 
0.0830 733 302 141 
0.0854 606 266 126 
6885 0.0745 Is Nas 273 128 
0.0818 868 281 Leal 
6886-2 0.0769 1,026 298 7a 
0.0830 582 215 124 


BPD UM|S HFUMNMN NM W WO WOM 


Sere) te) (02) LS) On C2) ony Sy Koy RY ee) ) 


PD OD WWE HP WHY WN OF 


WAN NO NUWD WO WA NY CF 


lRatio a/b refers to the ratio of the numbers indicated by 


a and b. 


*The SNR is computed from Oo, and 0_ based on averages with 
the outliers deleted. 
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are both accepted. Two of the seven records have maximum outlier ratios 
less than 5.3. This is 28 percent which is remarkably close to the 22 
percent derived above. 


The second measure of extremeness is the value of the SNR as based on 
o, and O_ computed without the outliers present. Other than noting 
that a three-sigma bound is often used in reliability to indicate unusual 
extremeness, no attempt will be made to interpret these results at this 
point in the research. 


One nonstatistical observation may be more significant than all the 
Statistical computations. This is the twofold fact that (a) the outliers 
always seem to fall in the maximum energy frequency range, and (b) the 
outlier lines are exceptionally separated from their neighbors, i.e., 
there is not a smooth transition with lots of small lines, some moderately 
large lines, and a few very large lines. It is more like two separate 
populations with no moderate range values between the two. This situation 
occurs in statistics in "gross error" or outlier questions. However, 
there is really insufficient statistical data to draw any firm conclusions. 


From an oceanographic viewpoint, the one or several outlier spectral 
lines might well dominate the waves present so that an aerial photo would 
show waves of that frequency proceeding in their particular direction. 
This is, of course, just conjecture since aerial photos for Hurricane 
Carla at that space-time location are not available. The other spectral 
lines present might be contributing noise and making the waves highly 
short-crested. 


The spectral line outliers may be some sort of resonant phenomena 
within the storm waves whereby energy tends to be concentrated on certain 
frequencies. Again, this is beyond the present research and is only 
noted in passing. 


In the situations where several outlier lines are present in the 
record, an investigation was made as to whether there could really be only 
one wave train present with the other lines showing up as leakage due to 
purely mathematical manipulations. Appendix B gives a derivation of the 
FFT leakage for a single cosine wave. It is shown that, at least in a 
gross sense, the leakage is delineated by the following formulas. Let 
the wave profile be given by: 


2T7mon 
Mea aces ( N - Ik (67) 


where Mo is not necessarily an integer. The approximate FFT spectral 
Finesse ate 20 <M +m << N, are given by the formula: 


2 sint(my - m) 2 
~ a“NAt { WG) } (68) 


T™(mpg - m) 
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From this, leakage would always be to the immediate neighboring lines. 
The outlier lines in Table 3 are always separated by one or more small 
lines in each case. Hence, one has to conclude that the cases with two 
or more outlier spectral lines cannot be explained by a single wave train 
with FFT leakage. 


XIII. WHY DOES CHI-SQUARED WORK FOR HURRICANE WAVES? 


One of the mysteries arising from the data is the surprisingly good 
probabilities arising from the chi-squared derivations. If linear wave 
theory was holding and the seas were Gaussian, this would be expected 
(Borgman, 1972, 1973). However, the hurricane waves were decidedly non- 
linear. The waves in many of the records have been plotted by computer 
and examined visually. The nonlinearities are really there. Why does the 
the chi-squared work so well? 


Investigation of this question led to a central limit theorem for 
dependent random variables which showed that the chi-squared relations 
hold exactly as N tends to infinity even for a non-Gaussian sea surface. 
The primary limitation is that water level elevations at the recorder 
separated by more than a certain constant time interval should be statis- 
tically independent of each other. The amount of the separation required 
to achieve independence is unimportant in the validity of the theorem 
although it will affect the speed of convergence to the asympototic 
result. Time sequences with above dependency properties are said to be 
"m-dependent ,"" in statistical terminology. 


A sufficient condition, then, for the central limit theorem to hold 
is that the probability density of the water level elevation measured 
from mean water level satisfies at least one of two "'tail'’ conditions. 
The density, f(m), should be such that there exists positive constants 
a, b, and c, and a positive integer n, such that: 


-b 
E)mcaalyy|inneg Wide Wbg| eas un: (69) 
or alternately, there exists a positive constant A _ such that: 

(iq) Ss OS Re sy] SA (70) 


The second condition is a special case of the first since if the second 
condition holds then c= A and any a, b, and n values will permit the 
first condition to be satisfied. 


Equation (69) is not an unreasonable restriction. For low seas, the 
sea surface has been found to be normally distributed. As the wave heights 
increase, a gamma density might be a reasonable guess as to the proper 
probability law. Both of these densities satisfy equation (69). In fact, 
from a practical viewpoint, no one seriously suggests that water level 
elevations can be infinite as required by the normal or the gamma densities. 
In fact, there will be a large value of A (e.g., A = water depth) such 
that equation (70) will hold. 
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The full details of the derivation for the m-dependent central limit 
theorem and the asymptotic chi-squared properties of the spectral esti- 
mates are given in Appendic C. One useful secondary result, summarized 
in item (C-8) of the appendix, is that the normed FFT coefficients, 


( Um f Vin 


Pubes AN | pee WMG) airs Lh Si oe ae (71) 
VNp,, At/2 VNp,, At/2 


asymptotically follow a multivariate probability law with zero mean and 
covariance matrix equal to an identity matrix. 


As an illustration of the asymptotic normality, 100 water level ele- 
vations were selected from the Hurricane Carla data number 6883. The 
elevations were taken 6 seconds apart from the beginning of the data 
until 100 elevations were obtained. Also, the 50 pairs of Fourier coef- 
ficients centered around the frequency associated with the largest spec- 
tral density value (Fig. 8) were tabulated from the computer listings. 
The water level elevations and the FFT coefficients were plotted on 
normal probability paper (after norming them) as shown in Figure 74. An 
examination of the figure shows that the FFT coefficients follow a rea- 
sonably straight line while the water level elevations exhibit a concave 
upward curve. The skewness of the water level elevations was computed to 
be 0.55 while the FFT coefficients only showed a skewness of 0.03. Thus, 
the FFT coefficients were reasonably close to normality even though the 
water level elevations were decidedly nonnormal. 


It would be interesting to explore the normalizing influence of the 
Fourier transform relative to the other wave records in Hurricane Carla. 
However, time did not permit that to be included in this investigation. 


XIV. IMPLICATIONS FOR DIRECTIONAL SPECTRUM RELIABILITY 


The asymptotic normality of the finite Fourier transform coefficients 
are extremely important relative to the estimation of the reliability of 
the directional spectrum and of other quantities such as energy diffrac- 
tion computed from the FFT coefficients. The FFT coefficients, at least 
for large data sets, can be taken as normally distributed and independent. 
These properties can be carried through the various formulas used in 
computing the particular quantity to obtain reliability measures for the 
quantity. 


XV. SUMMARY AND CONCLUSIONS 


1. The statistical properties of the spectral lines for 12 pieces of data 
measured during Hurricane Carla (8 to 10 September 1961) were examined 
in detail. The spectral lines were found to show negligible serial 
correlation. The spectral lines could be reduced to what appears to 
be random noise by subtracting from the lines the spectral density 
obtained by smoothing the lines with a moving average and then 
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dividing the positive and negative deviations by the local root-mean- 
square poSitive and negative deviations respectively. 


A simulation procedure based on the probability behavior of the random 
noise of the spectral line deviation was used to generate probability 
intervals for the spectral density estimates. 


The simulation probability intervals were found to agree reasonably 
well with that predicted from the chi-squared probability law. Al- 
though the chi-squared probabilities are expected to be valid for 
"low wave'' spectral densities, it was somewhat surprising to have 
them work for the nonlinear hurricane waves. 


An investigation of the above item led to the discovery that the 
finite Fourier transform coefficients will be approximately independ- 
ent and normally distributed for large data sets even though the water 
level elevations are not normally distributed. This face was derived 
from probability theory after assuming (a) m-dependence of the water 
level elevations about mean water level, and (b) a water level proba- 
bility density which is bounded for large y by the function 


aly|Me-bly| (a, b, and n are arbitrary constants). A central limit 
theorem was derived from these assumptions which led to the finite 
Fourier transform coefficients being asymptotically normally distrib- 
uted and independent as the data length tended to infinity. 


Under the conditions of the item above, the smoothed spectral esti- 
mates were found to be asymptotically chi-squared distributed pro- 
viding the true spectrum is constant over the smoothing interval. If 
the true spectrum is not constant, the spectral density estimates 
would still be reasonably close to having chi-squared behavior 
although there would be some deviations as evidenced by the deviations 
of the simulation results from the chi-squared results. 


The wave spectral lines exhibit members which appear unusually large 
as compared with the rest of the spectrum. The "outlier" lines 

always appear associated with the "peak energy" part of the spectral 
density. These unusual lines of energy might predominate the wave 
record so that a satellite picture would show waves of those frequen- 
cies riding on top of the rest of "noisy" wave combinations. However, 
this is just conjecture. 


The outlier spectral lines mentioned above do not persist for long 
periods of time. In the two cases where Hurricane Carla wave data 
were available 20 minutes apart, the outliers were found in one set 
of the data pair but not in the other. 


The conclusion of the study relative to reliability of directional 
spectrum estimation is that the finite transform coefficients may be 
taken as being serially uncorrelated and approximately normally 
distributed with zero mean, variance p,N At/2, and the appropriate 
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cross-variances. This asymptotic normality can be carried forward 
to provide reliability measures for quantities, such as the direc- 
tional spectrum, which are computed from the Fourier coefficients. 
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APPENDIX A 


EFFECTIVE WIDTH 


The effective width of a moving average smoothing function is defined 
to be the width of square pulse required to have the same area and same 
middle height. For example, consider the pulse: 


W(X) 


This pulse has area 


A -f w(x) dx 


and middle height w(0). The effective width would be the value of e.w. 
such that 


(e.w.) x w(0) =A 
or 


(e.w.) = A/w(0) 
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APPENDIX B 


EFFECT OF A DETERMINISTIC LINE ON THE SPECTRUM 
Suppose that i for O<n<N-1 is a cosine curve, 


i, = i eos [2m a | 5 (B-1) 


where mg is not necessarily an integer. What will be the effect on the 
spectral density? The fast Fourier transform of it gives: 


An = At y a cos [2m 7 ) 
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n=0 


Since, by the geometric series, 


DG rears Seely Sonia Nile d= Gs yWiG Ss se): 


it follows that: 


; i2m (mp-m) : ah 
hy = 25% [ete Le « gio 2 


e127 (mg+m) 
2 1 - et2™(mp-m)/N 97 12m (mp+m) /N 


(B-3) 
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ee aAt a oT (mg-m) et 7 (mp-m) eit (mp-m) 
n gi™(mg-m)/N =i (mp-m)/N _,im(mp-m)/N 


oi (my+m) iT (mg+m) _eit (mptm) 


eo LT (mg+m) /N ei ™ (mp+m) /N ti 7 tM (mgtm) /N 


aAt je my(1-q) - o} - ifm m1-4y} Sin MOmp-m) 
are sin ™(m)-m)/N 


e endl mp(1-<) - 6) - ifn ma-4y} Sin TMmp*m) 


sin ™(my+m)/N 


Soe sin T(mo-m) Ee sin 1(mp+m) 
. ant [er Byes 0 4 pCR) 0 (B-4) 
Sin ™(m-m)/N sin ™(mp+m)/N 
with 
d= TGs) 36 
= Tg (1-5 
1 
B = m™ (1-5) 5 (B-5) 
Hence, A, can be written in terms of real and imaginary parts as: 
aAt = T (mp-m) (a-B) Sin ™(mg+m) (a+B) 
= = | -———_ cos (a- + ————+——-_ cos (a+ 
An 2 \sriin 7 (mg-m) /N i sin ™(mg+m)/N 
. aAt fsin ™(mp-m) Sin 7 (ing+in) : 
+i = |= sin(a-8) - ——__" — _ ot : B-6 
2 E ™ (mj-m)/N one sin ™(m)+m)/N a 8)] Oe) 
It follows that: 
lan 2 r (Ger sin* 1(mg-m) i sin* T(mg+m) 
sin? 7 (mp-m) /N sin? 7 (mp+m) /N 
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> sin [7 (mp-m) | sin [1 (mg+m) | 


sin [ (mp-m) /N] sin [7 (mp+m) /N] {cos(a-B) cos (a+B) 


- sin(a-B) sin(a+8)}] 


ale 2 =. T (Mo-m) sin? T (mo+m) 
= — Ree oe 


sin? 7 (mp-m) /N sin? 7 (mo+m) /N 


2 sin [7 (mg-m) | sin [7 (mp+m) | 
ee eee 2 es 
sin [7 Gmg-m) /N] sin [7 (mp+m) /N] io «| ae 
4 N? sin? m(mp-m)/N  N* sin? t(mg+m)/N 


ae sin [m(mo-m)] sin [7 (mp+m) | 


1 
N2 sin [7 (mp-m) /N] sin [ (mp+m) /N] cos {2mmg (1-5) - 203 . (B-8) 


If mg is an integer with O<mp<N/2, 


ao NGA 5) alien Ss My 
P(fm) = (B-9) 


0 » otherwise. 


If N is much larger than either Mg-™ and Mg+m, which is usually 
the case in applications, then approximately, 


N sin 1(mg-m)/N ~ 1(mg-m) 
N sin T (my+m) /N x T (my+m) (B-10) 
1/N =~ 0 
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Hence, equation (B-8) can be written approximately as: 


B(E,) = ashe = T (™p-m) } ; 2 jee ied 


7 (mp-m) (m +m) 
Ags {= num {= noi) } cos (2 Tm - 26) . (B-11) 
FE m™ (mg-m) 7 (mp+m) 


If 20<mpgt+m << N, a reasonable assumption in many applications, then the 


first term inside the square brackets will dominate the other two. This 
fact is illustrated in the following computations, where 
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Hr are 
t= ee (mg ~ (B-12) 
1 7 (mg-m) 
in w(mg+m) ) 
sin T7(Mpyt+m 
— a {==} 2 HB) 2 0.0003 (B-13) 
2 T (mg+m) = LOO. ie 


fel ce ee onan 


T (mp-m) 7 (mp+m) 


il sin 7(mg-m) 
“70 7 T (mp -m) (B-14) 
My-M) = ORS Upper bound 
m qT, T, T, 
my 0.405 0.0003 0.020 
m,-1 0.045 0.0003 0.0068 
m,-2 0.016 0.0003 0.0041 
m,-3 0.0083 0.0003 0.0029 


9| 


Mpy-M, = OS 25 Upper bound 


m Th Dy TA 

m,+3 0.0067 0.0003 0.0026 
m +2 0.016 0.0003 0.0041 
n+l 0.090 0.0003 0.0096 
m a 0.0003 O.029 

m-1 0.032 0.0003 0.0057 
m,-2 0.010 0.0003 0.0032 
m -3 0.0048 0.0003 0.0022 


The main characteristics of equation (B-11) are determined by the first 
term and in a gross sense, 


2 
A _ aeNAt {== “tno-n) | i 
P(f,,) = mn Gian) (B-15) 


if 20<mp+m << N and |mp-m| << N. 
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APPENDIX C 


ASYMPTOTIC CHI-SQUARED PROPERTIES OF THE FFT SPECTRAL LINES 
FOR NON-GAUSSIAN, M-DEPENDENT WAVE TRAINS 


1. Basic Definitions and Assumptions. 


Let gs j = 0,41,+2,+3,... be water level elevations about mean 
water level. It will be assumed that {n,1 is a stationary second-order 


stochastic process which is not necessarily Gaussian and that the random 
sequence {njJ has the properties that, uniformly in n, 


E [nal = 0 (C-1) 


3 [ne SU (€=2) 


LOG Wee Maeee store vi S Onn A, Syeda siete | INES dis OG) is a finite 
sequence of the water level elevations starting at and terminating 
with NatN-] The time interval between water level elevation values is 
denoted by At. The finite Fourier transform coefficients are defined as: 


N-1 
a) = At » Y 


-i2T™mn/N 
nee 


(C-3) 


where i= v-l . The superscript N is attached to Am to indicate 
that the FFT coefficients are computed on the basis of a sequence of 
length N. 


It will be assumed that the sequence n. is m-dependent. That is, 


iy Tye ee eee Tyas? ang ifs 9 Merrersil? °° 


independent sets of random variables if b - c > m (Rosen, 1967). 


af no} are statistically 


The probability density for mn. will be assumed to satisfy either 
conditions (a) or (b) under item 5 given in the following. The FFT 
coefficients will be said to be degenerate if Dea 0 and hence 


Us = he 210: 


we) 


2. Motivation. 


Consider the set of nondegenerate coefficients: 


LOO 2) OD OD OD es GED GOD 
So eee Vingtl? Uingt2? Ving t2? 7 Mot Mo+T 


where Mo = goNS and <x» denotes the largest integer less than or 
equal to x. In the above, the constant a satisfies the inequality 
O<a<0O.5 and r is an integer constant with r>l. It will be shown that 
under the assumptions of item 1, the set S asymptotically is multi- 
variate normal as N+. If the true spectrum is constant over this set 
of r Fourier coefficients, then pP(fm)/p(fm) for the spectral lines 

for the frequencies spanned by the band will be asymptotically distributed 
as ol2 and the spectral density based on the average over the whole band 


will be asymptotically distributed as Ne |f Pies 
ZAG 


The first step, then, is to prove that the set S is asymptotically 
multivariate normal as N+. This requires the next two listed items. 


3. A Central Limit Theorem for m-Dependent Sums (Rosen, 1967). 
Consider the double sequence of random variables: 


(1) (1) (1) 
Xie haces Xe 


(2) (2) (2) 
Mg Mey booby Oe 


(N) (N) 
X5 Gsat0'o 8 xX 


(N) 
X 
l N 


> 


That is, the nth line of the array consists of a sequence of ky random 
variables. Let g (N be defined as: 


Ky 
eS Fa 
k=1 
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(i.e., the sum of the nth row). The central limit theorem is concerned 
with the conditions under which the probability law of 5 (N) converges to 


the normal probability law as N+. Other quantities which will be 
used in the theorem are: 
o2(s™) ) = Variance of g WN) (C-5) 
o2 = Variance of x ON) (C-6) 
k 
kN 
Ey d = probability density for ae 5 (Gey) 


Theorem: If the random variables in the same row of the array are 
m-dependent and if: 


(a) E(x) =0) VEonallyyke rand) IN); 


(b) aise SLM Nein Oe 
Ky 
(c) lim ewe eGo) ake S10) 4 soe GiGi, GSO. . 
ae kN 
kT ee 
and ky 
: 2 
(d) lim y ON <i Comm 
N?>o 
k=1 
then the probability law for 5 (N) converges to a normal probability law 


having zero mean and unit variance k tends to infinity. 


Comments: Rosen gives the theorem in a more general form by stating 
condition (c) in terms of the distribution.function and a Stieltjes inte- 
gral. However, the above form of condition (c) in terms of the probabil- 
ity density is sufficient for the present use. 


Proof: Given by Rosen (1967). 
4. Multivariate Central Limit Theorem (Rao, 1973). 
Let F, denote the joint distribution function of the k-dimensional 


random variable (a ae), Brana He Za mS Ueayoas eiaval JP the 


An 
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uP 'G 6.0 


distribution function of the linear function aT ae + yd 


n n 
+ Apa Also let F be the joint distribution function of a 
k-dimensional random variable a) z(2) ee 7 (kK) If for each vector 
elk > F,, the distribution function of i 71) + i z (2) algo SN z(k) 5 
An r 1 2 k 
then) Fa) aE. 


Proof: See Rosen (1967). 
5. Some Conditions for which (c) in the Theorem of Item 3 Holds. 


Suppose that c are constants uniformly bounded in n and k, 


nk 


eS ee il 


and the probability density of Ue is denoted by ee Og) Tifa: 


k 


(a) there exists positive constants a, b, and c such that 
=D) gy 
> 
eu Opus ee 2) en llbalene 


uniformly in n and k, or 


(b) there exists a positive constant A_ such that 
el) SO 5 ae yea 


uniformly in n and k. 


Then for any e€>0, 


RS 2 2 pas 
L = lim if xe £2 dx = 0 
k= ed | ele 


Proof: The probability density for Mo expressed in terms of £2) is: 


$0) = ( lene? le (ene 7 7 | 
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Hence, L may be written in terms of y integration as: 


n 
L = lim c 2 |c | fs “nk d 
nk” |°nk| *nk Tx 
DAY k=1 Ca 5S¥ n z 
nk Yn 
—|>e 
vn 
n 
— hoe ik 2 2 
= ee = », cat | Pon OGY 
k=1 ly[peyn/[c., | 


Let B be the uniform bound on {e_,} : 


(a) Proof of the conclusion under hypothesis (a) above: for a>0, 
and fixed n,a, and b, define G(a) as: 


foo} foe) 


G(a) -f yeaynen olay, = — i oe ae 
b / 
a ba 


Clearly G(a) is a monotone decreasing function of a and G(a)*0 as 
aro, 


Hence, including both tails of eewy) 5 


a n 
1 2 
n » Ck a (O) chy S = \: ee nee bY ay 
evn 0 = eal 
k=1 Palate i k=1 yo 
nk [cx 


as n tends to infinity. 


(b) Proof of the conclusion under hypothesis (b) above: BY) 
Satisfies hypothesis (a) with C =A, n= 0, a = 1, and b = 1. 


SI, 


6. <A Useful Lemma. 


etn iee - C5 > Cyy-ss be a sequence of constants and Wi> W,, Wz, 
be a sequence of random variables. Suppose that 


lim 
Co Ses 
I: rene Tip) 


and that there exists a random variable Z with finite variance such that 
the probability law for c,W, converges to the probability law for Z as 
n-* © . Then the probability law for cW, also converges to that for Z 


als) IM =e 2" g 


Proof: The proof is a straightforward application of a relation given by 
Rao (1973). In his terminology, the above lemma may be stated as: 


L L 
CpW,> Z implies cW,> Z 


His relation states that this holds provided, for any e€ > 0, 


iam P(|c W. - cW [>] = 0 
nn n 


no7o 


P 
(In Rao's notation this would be stated as |coaWn = cW, | > 0). Now 


See W W 
F(lenk, = Mal > | = &|[m,) > =| < a 
nn n n [en-¢| 22 


by the Tchebichev inequality (Loeve, 1960). Since Var(W,) > Var(Z) 
and (c, - c)2 >0 as n+o© , it follows that: 


lim 
aN P| lentin = cWnll > | = 0 
as required. Hence, the lemma is proven. 


7. Asymptotic Normality of Linear Combinations of Nondegenerate FFT 
Coefficients. 


Let Ajo tor s = 1,2 sand) mp + 1 semis mp + © betany sequencesor 


bounded constants and define: 


roy, y ON) POnt y ON) 
(N) _ m \ m 
ui iz » AS ‘ » m2 
Ty NPAGp a2 = yNeepeyi2 
m=Mo+1 m=Mg+1 
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where the various terms are defined in item 1, and the assumptions listed 
there hold. Let o2(T(N)) be the variance of T(N) , Then, 


My+r 


ee co eu) ale ce i 2, ) 


m=Mg+1 


The probability law of 


converges to a normal probability law with zero mean and unit variance as 
N tends to infinity. 


Proof: The terms wee and ve are asymptotically uncorrelated with 


each other and with other pairs (yu 5 vod ) provided 0 <m< N/2 
and 0 < m' < N/2 (Borgman, 1973). The asymptotic variance of both 


mee and Vinge is PN At/2 (Borgman, 1973). Hence, U,/ WN Atp_/2 and 


vA Atp,/2 have unit variance. It follows that: 


Motr 


san) 0 A ee 


m=M +1 


as N > o (Freund, 1971). 


After substitution for Un and V in terms of Y_, 7) can be 
expressed as: ui u 


i k-1 : k- 
7 ON) : 5 mop 5 28 oa cos ann) +d, sin eimkol)) 
o(7)) k=l YN o(1)) m=mg+1 yp, At/2 
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Let 


mtr (a ee 27m (k-1) Be eat aro) ) 
5 ml m2 N 
eae 
Nk N 
m=m+1 o(r' )) [> At/2 
N 
ne Snk vas ) 
ee qh way 
and 
ky = N 
Then, 
7 WN) 


1:0) ) Dm 


is in the form required for the application of the theorem in item 3. 


Since Eve) = Elna] = 0 , condition (a) is satisfied. The 
division by o(T(N)) satisfied condition (b). The variance of X(N) 


is given by: 


kN N 
But , 
Myjt+r 
0 
| | Z [nal i nal im 
“nkl = », (N) een 
m=mp*l o(t ) /P.. At/2 
so, 
N 5 var(y.)) N (N) 
2 —_—S 2 = 2 
», OrN 5 Ser var( yf ) B 
k=11 k=1 
which is bounded as N>*©. This verifies condition (d). 
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Finally condition (c) will be satisfied if the probability density 
for the water level elevations obey the bounding conditions (a) and (b) 


in item 5. This was assumed in item 1. (Note: B= lim By + 1 will 
> © 
provide the uniform bound for cok needed in item Sal 


Thus, all conditions are satisfied and POND 75 cp ON) converges in 
law to a zero mean, unit-variance normal probability law. Since 


Mp)tr 
2 2 A \ ae 
PE) Ces) 0 
m=M+1 


as N>+o , it follows from item 6 that the probability law of TN) jg 
converges to zero mean, unit-variance normal probability law. This 
completes the proof. 


8. Asymptotic Normality of the FFT Coefficients. 


Let (ue ‘ vo) for m= Mot, Mot2, slots INO +e De asiet ot: 


nondegenerate FFT coefficients. Then the multivariate probability law 
for 


y™ y) 
m 


( Beech Van eee ae 
Jwp, at/2 Np, t/2 


) pe casa 


will be a multivariate normal with covariance matrix equal to an identity 
matrix and mean vector having all components equal to zero. 


Proof: In items 4 and 7. 


9. Asymptotic Chi-Squared Distribution for Spectral Estimates. 


For the range of m-values, mp+l < m < mp+r, suppose that Pease 


is a nonzero constant. Let: 


+ 


p=s y (u2 + v2) /n ae 


m=Mg+1 
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Then, 2rp/p is asymptotically a chi-squared random variable with 2r 
degrees of freedom. 


Proof: 
My*r U 2 V 2 
rir Spoon} (igen not rears 
At At 
mM=Mg+1 
or 
Mg*r U 2 V 2 
P oy LS YNP At/2 WNp At/2 
M=Mo 


By item 8, the terms squared in the sum are asymptotically independent, 
zero mean, unit-variance normal random variables. Hence, the right-hand 


side is asymptotically chi-squared with 2r degrees of freedom (Freund, 
1971). 
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